In [1]:
import dicom
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
In [2]:
dcm_plan = dicom.read_file("RP1.3.6.1.4.1.2452.6.728760501.1211531093.3851412398.2396712583.dcm", force=True)
# dcm_plan
In [3]:
dcm_plan.ApplicationSetupSequence[0].dir()
Out[3]:
['ApplicationSetupNumber',
 'ApplicationSetupType',
 'ChannelSequence',
 'TotalReferenceAirKerma']
In [4]:
dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0].dir()
Out[4]:
['BrachyControlPointSequence',
 'ChannelLength',
 'ChannelNumber',
 'ChannelTotalTime',
 'FinalCumulativeTimeWeight',
 'NumberOfControlPoints',
 'RefdROINumber',
 'RefdSourceNumber',
 'ReferencedROINumber',
 'ReferencedSourceNumber',
 'SourceApplicatorID',
 'SourceApplicatorLength',
 'SourceApplicatorName',
 'SourceApplicatorNumber',
 'SourceApplicatorStepSize',
 'SourceApplicatorType',
 'SourceMovementType',
 'TransferTubeLength',
 'TransferTubeNumber']
In [5]:
dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0].BrachyControlPointSequence[0]
Out[5]:
(300a, 0112) Control Point Index                 IS: '0'
(300a, 02d2) Control Point Relative Position     DS: '30.0000000000000'
(300a, 02d4) Control Point 3D Position           DS: ['-36.132858', '11.474540', '75.310880']
(300a, 02d6) Cumulative Time Weight              DS: '0.00000000000000'
In [6]:
dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0].BrachyControlPointSequence[0].dir()
Out[6]:
['CP3DPosition',
 'CPIndex',
 'CPRelativePosition',
 'ControlPoint3DPosition',
 'ControlPointIndex',
 'ControlPointRelativePosition',
 'CumulativeTimeWeight']
In [7]:
dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0].BrachyControlPointSequence[0].ControlPoint3DPosition
Out[7]:
['-36.132858', '11.474540', '75.310880']
In [8]:
dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0].BrachyControlPointSequence[0].CumulativeTimeWeight
Out[8]:
'0.00000000000000'
In [9]:
number_of_channels = len(dcm_plan.ApplicationSetupSequence[0].ChannelSequence)
number_of_channels
Out[9]:
18
In [10]:
number_of_dwells_in_channel = len(
    dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0].BrachyControlPointSequence)
number_of_dwells_in_channel
Out[10]:
10
In [11]:
all_dwell_positions_by_channel = []
all_dwell_positions_flat = []

for i in range(number_of_channels):
    BrachyControlPointSequence = (
        dcm_plan.ApplicationSetupSequence[0].ChannelSequence[i].BrachyControlPointSequence)
    number_of_dwells_in_channel = len(BrachyControlPointSequence)
    
    dwells_positions = []
    for j in range(number_of_dwells_in_channel):
        if len(BrachyControlPointSequence[j].dir("ControlPoint3DPosition")) != 0:
            dwells_positions.append(
                BrachyControlPointSequence[j].ControlPoint3DPosition)
    
    all_dwell_positions_by_channel.append(dwells_positions)
    all_dwell_positions_flat += dwells_positions
    
all_dwell_positions_by_channel
Out[11]:
[[['-36.132858', '11.474540', '75.310880'],
  ['-36.132858', '11.474540', '75.310880'],
  ['-36.132858', '11.474540', '65.315258'],
  ['-36.132858', '11.474540', '65.315258'],
  ['-36.237490', '11.474540', '55.317446'],
  ['-36.237490', '11.474540', '55.317446'],
  ['-36.237490', '11.474540', '45.321398'],
  ['-36.237490', '11.474540', '45.321398'],
  ['-36.215044', '11.265048', '35.336656'],
  ['-36.215044', '11.265048', '35.336656']],
 [],
 [['3.758867', '11.228837', '75.297996'],
  ['3.758867', '11.228837', '75.297996'],
  ['3.759700', '11.230399', '65.298029'],
  ['3.759700', '11.230399', '65.298029'],
  ['3.717219', '11.310051', '55.303095'],
  ['3.717219', '11.310051', '55.303095'],
  ['3.716297', '11.311779', '35.303131'],
  ['3.716297', '11.311779', '35.303131']],
 [],
 [['43.944967', '11.474540', '65.256192'],
  ['43.944967', '11.474540', '65.256192'],
  ['43.945226', '11.474540', '55.260574'],
  ['43.945226', '11.474540', '55.260574'],
  ['44.049858', '11.474540', '45.262763'],
  ['44.049858', '11.474540', '45.262763'],
  ['44.189367', '11.474540', '35.266653'],
  ['44.189367', '11.474540', '35.266653']],
 [],
 [],
 [],
 [],
 [],
 [['-16.357477', '-9.423887', '87.754693'],
  ['-16.357477', '-9.423887', '87.754693'],
  ['-16.357477', '-9.423887', '77.754693'],
  ['-16.357477', '-9.423887', '77.754693'],
  ['-16.357477', '-9.423887', '67.754693'],
  ['-16.357477', '-9.423887', '67.754693'],
  ['-16.357477', '-9.521285', '47.756594'],
  ['-16.357477', '-9.521285', '47.756594']],
 [],
 [['24.072189', '-9.277403', '77.770644'],
  ['24.072189', '-9.277403', '77.770644'],
  ['24.072189', '-9.277403', '57.770644'],
  ['24.072189', '-9.277403', '57.770644'],
  ['24.071354', '-9.277635', '37.776129'],
  ['24.071354', '-9.277635', '37.776129']],
 [],
 [],
 [],
 [],
 [['38.818276', '-19.775444', '87.761693'],
  ['38.818276', '-19.775444', '87.761693'],
  ['39.143063', '-19.721680', '67.772517'],
  ['39.143063', '-19.721680', '67.772517'],
  ['39.252303', '-19.612684', '57.783277'],
  ['39.252303', '-19.612684', '57.783277'],
  ['39.305804', '-19.676885', '47.784693'],
  ['39.305804', '-19.676885', '47.784693'],
  ['39.359840', '-19.613848', '37.794699'],
  ['39.359840', '-19.613848', '37.794699']]]
In [12]:
dwell_positions = np.array(all_dwell_positions_flat).astype(float)
dwell_positions
Out[12]:
array([[-36.132858,  11.47454 ,  75.31088 ],
       [-36.132858,  11.47454 ,  75.31088 ],
       [-36.132858,  11.47454 ,  65.315258],
       [-36.132858,  11.47454 ,  65.315258],
       [-36.23749 ,  11.47454 ,  55.317446],
       [-36.23749 ,  11.47454 ,  55.317446],
       [-36.23749 ,  11.47454 ,  45.321398],
       [-36.23749 ,  11.47454 ,  45.321398],
       [-36.215044,  11.265048,  35.336656],
       [-36.215044,  11.265048,  35.336656],
       [  3.758867,  11.228837,  75.297996],
       [  3.758867,  11.228837,  75.297996],
       [  3.7597  ,  11.230399,  65.298029],
       [  3.7597  ,  11.230399,  65.298029],
       [  3.717219,  11.310051,  55.303095],
       [  3.717219,  11.310051,  55.303095],
       [  3.716297,  11.311779,  35.303131],
       [  3.716297,  11.311779,  35.303131],
       [ 43.944967,  11.47454 ,  65.256192],
       [ 43.944967,  11.47454 ,  65.256192],
       [ 43.945226,  11.47454 ,  55.260574],
       [ 43.945226,  11.47454 ,  55.260574],
       [ 44.049858,  11.47454 ,  45.262763],
       [ 44.049858,  11.47454 ,  45.262763],
       [ 44.189367,  11.47454 ,  35.266653],
       [ 44.189367,  11.47454 ,  35.266653],
       [-16.357477,  -9.423887,  87.754693],
       [-16.357477,  -9.423887,  87.754693],
       [-16.357477,  -9.423887,  77.754693],
       [-16.357477,  -9.423887,  77.754693],
       [-16.357477,  -9.423887,  67.754693],
       [-16.357477,  -9.423887,  67.754693],
       [-16.357477,  -9.521285,  47.756594],
       [-16.357477,  -9.521285,  47.756594],
       [ 24.072189,  -9.277403,  77.770644],
       [ 24.072189,  -9.277403,  77.770644],
       [ 24.072189,  -9.277403,  57.770644],
       [ 24.072189,  -9.277403,  57.770644],
       [ 24.071354,  -9.277635,  37.776129],
       [ 24.071354,  -9.277635,  37.776129],
       [ 38.818276, -19.775444,  87.761693],
       [ 38.818276, -19.775444,  87.761693],
       [ 39.143063, -19.72168 ,  67.772517],
       [ 39.143063, -19.72168 ,  67.772517],
       [ 39.252303, -19.612684,  57.783277],
       [ 39.252303, -19.612684,  57.783277],
       [ 39.305804, -19.676885,  47.784693],
       [ 39.305804, -19.676885,  47.784693],
       [ 39.35984 , -19.613848,  37.794699],
       [ 39.35984 , -19.613848,  37.794699]])
In [13]:
from plotly.graph_objs import *
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
init_notebook_mode()

dcm_struct = dicom.read_file("RS1.3.6.1.4.1.2452.6.209502505.1199599685.3559183269.2231927106.dcm", force=True)

structure_names = [item.ROIName for item in dcm_struct.StructureSetROISequence]
structure_names
Out[13]:
['Prostate PTV',
 'Urethra',
 'Rectum',
 'Seed Apex',
 'Seed Base',
 'Seed Mid',
 'Avoidance1',
 'Avoidance2',
 'Solid Water',
 'Phantom',
 'Detector']
In [14]:
def pull_structure(number):
    structure_names = [
        item.ROIName for item in dcm_struct.StructureSetROISequence]
    
    contours_by_slice_raw = [
        item.ContourData for item in dcm_struct.ROIContourSequence[number].ContourSequence]
    x = [np.array(item[0::3]) for item in contours_by_slice_raw]
    y = [np.array(item[1::3]) for item in contours_by_slice_raw]
    z = [np.array(item[2::3]) for item in contours_by_slice_raw]
    
    print("Loaded {}".format(structure_names[number]))
    return x, y, z


def display_structures_with_dwells(list_of_structures, dwells, colour_list=None):
    dicom_structure_names = np.array(
        [item.ROIName for item in dcm_struct.StructureSetROISequence])
    combined_trace = []
    
    for i, structure in enumerate(list_of_structures):
        if colour_list is None:
            colour = 'black'
        else:
            colour = colour_list[i]
        
        reference = (dicom_structure_names == structure)
        if np.all(reference == False):
            raise Exception("Structure not found (case sensitive)")

        index = int(np.where(reference)[0])
        x, y, z = pull_structure(index)
        
        for i in range(len(x)):
            trace = Scatter3d(
                x=np.append(x[i], x[i][0]), 
                y=np.append(y[i], y[i][0]), 
                z=np.append(z[i], z[i][0]),  
                mode='lines', line=Line(color=colour, width=3))

            combined_trace.append(trace)

    trace = Scatter3d(
        x=dwells[:,0], 
        y=dwells[:,1], 
        z=dwells[:,2], mode='markers')
    combined_trace.append(trace)
            
    iplot(Figure(
            data=Data(combined_trace),
            layout=Layout(
                showlegend=False,
                width=900, height=800
            )
        ))

Final diagram

In [15]:
display_structures_with_dwells(
    ['Prostate PTV', 'Urethra', 'Rectum'], dwell_positions,
    ['Red', 'Blue', 'Green']
)

Loaded Prostate PTV
Loaded Urethra
Loaded Rectum
In [ ]: